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ABSTRACT 



Context. Recent observations from NASA's Kepler mission detected the first planets in circumbinary orbits. The question we try to 
answer is where these planets formed in the circumbinary disk and how far inside they migrated to reach their present location. 
Aims. We investigate the first and more delicate phase of planet formation when planetesimals accumulate to form planetary embryos. 
Methods. We use the hydrodynamical code FARGO to study the evolution of the disk and of a test population of planetesimals 
embedded in it. With this hybrid hydrodynamical-N-body code we can properly account for the gas drag force on the planetesimals 
and for the gravitational force of the disk on them. 

Results. The numerical simulations show that the gravity of the eccentric disk on the planetesimal swarm excites their eccentricities 
to values much larger than those induced by the binary perturbations only within 10 AU from the stars. Moreover, the disk gravity 
prevents a full alignment of the planetesimal pericenters. Both these effects lead to large impact velocities, beyond the critical value 
for erosion. 

Conclusions. Planetesimals accumulation in circumbinary disks appears to be prevented close to the stellar pair by the gravitational 
perturbations of the circumbinary disk. The observed planets possibly formed in the outer regions of the disk and then migrated inside 
by tidal interaction with the disk. 



Key words. Planets and satellites: dynamical evolution and stability — Planet-Disk interactions - 



Methods: numerical 



1. Introduction 

Most Sun-like stars are believed to form as gravitationally bound 
pairs. Estimating the prevalence of planets in binary systems is 
then a relevant issue when computing the fraction of stars with 
planets. Planets in circumstellar orbits (S-type orbits) around 
one of the components of a binary system have been found 
in both close (~20AU) and wide binary systems confirming 
that the binary perturbations may not be able to prevent planet 
formation even if they can significantly affect the course of 
it. This last issue has been explored in detail using numeri- 
cal simulations, most of them modeling the intermediate stage 
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2008; Ma rzari et al.M2.012b . These studies have shown that this 
phase might be the most sensitive to binarity effects, because 
mutual impact velocities can be increased to values that may 
threaten the formation of large bodies. Crucial in this phase is 
the evolution of both the eccentricity and perihelion longitude 
of the planetesimals under the coupled action of the compan- 
ion gravity and gas drag force. A size-dependent phasing of the 
orbits develops over the timescale of the secular perturbations 
of the binary, leading to high (and accretion inhibiting) colli- 
sion velocities a mong any planetesima l population with even a 
small size spread (The bault et al.l l2006). Mechanisms that could 
come to the rescue of the planetesimal accumulation process are 
a small inclination between the circumstellar disc and the bi- 



nary plane dXie et al ] 1201 Obi) , or an outward migration of pro- 
toplanetary embry os formed in safer regions closer to the star 
dPavne et al.ll2009l) . or planetesimal growth through the sweep- 
ing of small collisiona l fragments dPaardekooper & Leinhardtl 
12010b IXie et al.ll2010al) . or the fact that the binary was initially 
wider and was compacted by stellar encounters d uring the early 
evoluti on of the stellar cluster it was born in dThebault et al.l 
l2009bl) . Another, more radical solution would be that planet for- 
mation proceeds through a different channel in close binaries, 
a hypothesis that could be supported by the fact the exoplan- 
ets found in close bi naries have diff erent properties than those 
around single stars (lDuchendl20To|) . For a more detailed dis- 
cussion on circu mprimary planet formation in binaries, see the 
recent review by iThebauld (1201 lb . 

Recently, the KEPLER team announced the discovery of an- 
other category of exoplanets in binaries, i.e, transiting circumbi- 
nary planets, which are planets moving on P-type orbits circum- 
venting both the st ars. The first of su ch planets to be discovered 
was Kepler- 16 b dDovle et al.l|2011l) followed by Kepler-34 b, 
Kepler-35 b (IWelsh et al.ll2012l) and Kepler-47 b,c the first cir- 
cumbinary multi-planet system (lOrosz et al 1 120121) . The forma- 
tion of planets in circumbinary P-type orbits is very different 
compared to that in S-type orbit (around a component of the 
pair) since the secular perturbations of the companion star on a 
planetesimal swarm have different intensity and form compared 
to that produced by an external perturber. As a consequence, the 
outcomes of numerical modeling of planetesimal evolution in 
S-type orbits cannot be applied to study planet formation in cir- 
cumbinary orbits. 
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At present, two numerical studies have been performed 
to estimate where in the initial circumbinary disk around 
Kepler-16 planetesimal accumulation can proceed, even if per- 
turbed and possibly at a slower pace, towar ds the accumula- 
tion of large planetary em bryos. Both studies (lMeschiarill20 1 2t 
iPaardekooper et ai1l2012l) used an N-body approach to compute 
the trajectories of a large number of planetesimals perturbed by 
the non-spherically symmetric gravit ational field of the stellar 
pair. The first paper (lMeschiarill20 1 2b focused on the long term 
effects of the secular perturbations finding potential accretion- 
friendly zones within 1.75 AU from the star pair and beyond 4 
AU. In these regions the mutual impact velocity between plan- 
etesimals were, in the majority of collisions, lower than the 
threshold velocity causing disruption. Their eccentricity distri- 
bution is centered around the the forced component of the secu- 
lar perturbations of the binary given by 



e f = -(1 - 2ju) — e B , 



(1) 



where // = molimx + m.2) is the binary mass ratio, a# and eg 
the semimajor axis and eccentricity of th e binary system and a is 
the se mimajor axis of the planetesimal dMoriwaki & Nakagawal 
12001 . 

In a subsequent paper, IPaardekooper et aD (1201 2l) re- 
analysed planetesimal accumulation around Kepler-16 and also 
explored two other systems, Kepler-34 and Kepler-35. They fo- 
cused on the effects of short term eccentricity perturbations on 
the planetesimal motion. They also considered the possible reac- 
cumulation of small fragments, produced by the shattering of 
planetesimals impacting at high relative velocity, onto the largest 
intact planetesimals. Even if th is second-generation accretion 
helps. IPaardekooper et al.l (120 1 2b conclude that planet formation 
is not possible at the present location of the planet, but may 
however be effective beyond twice the present semi major axis 
(a - 0.7 AU). Both the previous ly mentioned studies dMeschiaril 
l2012t IPaardekooper et aDl2012l) claim that, in spite of the com- 
bined perturbations of the binary star system and of the friction 
from the gas of the circumbinary disk, planet formation seems 
to be possible beyond a distance of about 2-4 AU from the bari- 
center of the pair, which is further out than the present location 
of the exoplanets but still relatively close to the binary. 

The basic assumption of an N-body approach to the problem 
of planetesimal accumulation requires that the gas disk in which 
the planetesimals are embedded is axisymmetric. This approxi- 
mation may not be a good one when dealing with circumbinary 
disks since the gravitational perturbations of the sta r pair may 
strongly affect the disk shape. As already shown in (|Klev et al] 
teOOalMarzari et al.l2009bll2012tlMuller & Klevl2012E the disk 
may become eccentric and be perturbed by strong spiral waves. 
In that paper, the parameters of the system were similar to those 
of Kepler-16 but the disk was truncated at 3 AU and the isother- 
mal approximation was used. In this paper we adopt a more com- 
plete approach modeling the circumbinary disk that gave origin 
to the planet in Kepler-16 as a radiative disk extending out to 10 
AU from the baricenter of the two stars. The radiative model is 
better suited to describe the earlier stages of the disk evolution 
when it is massive and probably optically thick. It has also been 
shown (Marza ri et al ] l2012t iMuller & Klevll20T2h that radiative 
disks, when perturbed, develop overall shapes and internal struc- 
ture that differ from those of locally isothermal disks, in particu- 
lar concerning the disk eccentricity and the propagation of spiral 
waves which may significantly perturb planetesimal trajectories. 
For this reason, in order to obtain an accurate modeling of the 



disk structure, in our simulations we solve an energy equation 
that includes the viscous heating of the disk and radiative losses. 
On the basis on the above mentioned references we expect that 
a detailed treatment of the disk thermodynamics is more rele- 
vant for the evolution of the disk eccentricity than the inclusion 
of the effects of self-gravity. In addition, considering a larger 
disk allows a better handling of the gravitational perturbations 
of the disk on the planetesimal orbits. If the disk is more mas- 
sive, since we model a larger portion of it, the perturbations due 
to a potentially uneven mass distribution, due to the building up 
of an eccentric shape, are considerably stronger. 

In this paper we compute the trajectories of a large number 
of planetesimals perturbed by the primordial circumbinary disk. 
Our simulations focus on the Kepler-16 system and show that 
the gravitational perturbations of the disk excite large eccentric- 
ity values in the pla netesimal swarm w i th a m echanism similar 
to that described by iNelson & Gressell (120101) . These large ec- 
centricities, up to 0.4 and beyond, are able to prevent the onset 
of accumulation within 10 AU from the stars and maybe also 
farther out. The present planet observed in the system possibly 
formed in t he outer regions of the disk an d migrated inside as 
described in lPierens & Nelson! (I2007U2008I) . The migration pos- 
sibly occurred in all the circumbinary systems mentioned above 
and is strongly suggested by the coincidence between the semi- 
major axis of the planets and the location of the internal stab ility 
limit which can be derived from lHolman & Wiegertl (11999b . As 
a consequence, our results do not suggest that planet formation 
is not possible, but that it had to occur in the outer regions of the 
disk where the perturbations of the star pair are less effective in 
producing non-axisymmetric density perturbations on the disk. 

In Section 2 we describe the numerical model and in Section 
3 we describe out results. Section 4 is devoted to the discussion 
of the results. 

2. The hybrid algorithm modeling the evolution of 
the disk and planetesimals 

To compute the trajectories of planetesimals and, at the same 
time, the evolution of the gaseous disk under the perturbations 
of the bin ary, we have u sed the two-dimensional numerical code 
FARGO dMassetll2000l) modified to fit the problem. The hydro- 
dynamical equations are solved in a cylindrical coordinate sys- 
tem centered on the baricenter of the binary. We focus on the 
Kepler-16 system where the mass of the primary is Mi = 0.69 
M , that of the secondary Mi = 0.20 M G , the binary semima- 
jor axis as = 0.224 AU and the eccentricity eg = 0.159. The 
two stars are evolved on a fixed Keplerian orbit neglecting any 
change in the binary system due to gas accretion by the stars and 
any momentum exchange with the disk. T his choice is justified 
by the models of lPierens & Nelson! J2007h suggesting that a sys- 
tem made of the binary and the disk reach a near-stationary state 
after some evolution. Incidentally, these authors used binary pa- 
rameters actually close to those derived for Kepler-16. 

The grid used in our calculations to model the disk has 
N r = 256 radial zones and A^ s =512 azimuthal zones, and an 
arithmetic spacing is used along the radial direction i.e. the ra- 
dial distance is divided in equal size intervals. All the simula- 
tions are carried out inclu ding an energy equ ation of the form 
dBaruteau & MassetEool iMarzari et al.ll2012l) 

_ + V • (ev) = - P V ■ v + (2+ sc - Q- cml + AeV 2 log(p/57) (2) 

where e = p/(y - 1) is the thermal energy density, y = 1.4 is 
the adiabatic index, and v denotes the gas velocity. In the equa- 
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tion we do not include the effects of stellar irradiation. The term 
<2 + isc is the heating term due to the viscous heating while the 

cooling term 2c 00 i ' s assume d to be 2<x SB 7' 4 ff , where cr SB is the 
Stefan-Boltz mann constant and T e g is the effective temperature 
estimated as (lHubenvlll990h 

T 4 eff = r 4 /r eff , 

for an effective optical depth 



(3) 



Teff 



3r V3 



1 

4t' 



(4) 



The vertical optical depth, r, is approximated as r = kZ/2, 
wher e for the Ro s seland mean opac ity, k, we adopt the formu- 
lae in lBell&Linl d!994l) . Following IPaardekooperet al l (1201 ll) , 
we also model thermal diffusion as the diffusion of the gas en- 
tropy, s, defined as s = ft.(y - log(p/E y ). This corresponds 
to the last term in the right-hand side of Eq. (f2}, where A is a 
constant thermal diffusion coefficient. Throughout this study, we 
adopt A = 10 -6 in code units. The initial aspect ratio h = H/r is 
constant all over the disk and is set to 0.05. A constant shear 
kinematic viscosity of 10~ 5 (normalized units), which corre- 
sponds at about 5 AU within t he dis k to an a value of about 
2.5 x 10 3 (S hakura & Sunvaevld 19731) ). is used and open bound- 
ary conditions are adopted with standard outflow at both the in- 
ner and outer edge. The initial temperature profile of the disk 
is computed as T(r) = T$r~ l where the value at 1 AU is set 
to Tq = 630 K. This value is der ived following the approach 
described in iMarzari et al.l (1201 2l) and it depends on the ini- 
tial choice of the aspect ratio h. The Toomre-Q parameter is 
quite large in the inner disk parts, where the binary's pertur- 
bation is the strongest. Its radial dependence is approximately 
Q(R) ~ 100 x (R/lAUy 3/2 and a value of about 10 is measured 
at R = 5 AU. 

We neglect in our model the apsidal precess ion of the binary 
due to its interaction with the disk. According to Rafi kovl(l2012l) . 
the period of the binary precession due to an axisymmetric disk 
is given by 



T g , = *ir\ 



(M B \ 


( 1/2 5/2 \ 


\M d j 


a B 3 <pn B/ 



(5) 



where is the disk mass, Mb the summ of the star masses, 
n B the mean motion of the binary, r,„ the inner border of the disk 
and r the outer one, and is a constant whose value depends 
on the ratio of a B and r,„ and on the mass ratio and it can be ap- 
proximated by 0.5. Using this equation we find that T a ~ 3 x 10 4 
yrs (after some initial evolution, the mass of the disk settles to 
Md = 2.75 x 1O~ 2 M ). This is longer than the timescale over 
which the planetesimal eccentricity grows i.e. 10 3 yrs. We ex- 
pect that the binary precession might have an effect on the long 
term evolution of the system but not on the short timespan we 
are covering with the model. In addition, since the disk in our 
simulations is eccentric, Eq. Q may not be very precise being 
derived under the assumption of an axisymmetric disk. As a con- 
sequence, a full numerical approach is needed s ince the formul a 
for the potential of a axisymmetric disk, used in iRafikovl (|20 1 2|) . 
cannot be applied. We plan to explore this effect in the future 
hoping in an increase in computing power. 

The forces acting on the planetesimals include the binary 
gravitational force, the gravitational force exerted by the disk 
and the g as drag force. The latter is calculated in the Stokes 
regime as dAdachi et al.ll 1976b 



where v is the relative velocity between the planetesimal and 
the gas and the parameter k (Kary et al. 1993) is given by 



k = 



3p g C D 
8p„Rp 



(7) 



where Co is the aerodynamic drag coefficient for objects 
with large Reynold's number, p p andR p are planetesimal density 
and radius, respectively, and p g is the gas density in the midplane 
derived from the surface mass density S(r) through the relation 



Pg = Z(r)/[(2n) 1/2 H] 



(8) 



F = kvv 



(6) 



with H the disk scal e height dGiinther & Kley|2002l) . Respect 
to (IMarzari et al 1 120081) we adopt a more extended disk ranging 
from 0.5 to 10 AU from the binary system with a superficial 
density profile S = L^r ~ 1 ^ 2 where So is the density at 1 AU set 
to So = 2.5 x 10~ 4 in normalized units, i.e. ~ 2.2 x 10 3 g/cm 2 , 
compatible with Minimum Mass Solar Nebula density. 

To refine the computation of the drag force acting on the 
planetesimals we extrapolate the density and velocity (p, v) of 
the gas at the planetesimal location with a bilinear fit from the 
values at the bord ers of each grid cell. T his is a refinement of the 
algorithm used in lMarzari et al.l(l2008l) . 

3. The results 

In this section we compute the orbital evolution of the planetes- 
imals and of the disk. Before including the planetesimals in the 
model, we evolve the disk for 10 5 binary revolutions (approx- 
imately 10 4 yr). We then restart the simulation including 400 
planetesimals on initially circular orbits with semimajor axis 
equally spaced from 1 to 8.8 AU and we compute their orbits 
and the disk evolution for 10 4 yrs. 



3. 1 . The disk shape 

In the upper panel of Fig. [2] we show isocurves of the surface 
density distribution of the disk confirming its eccentric shape 
responsible for strong gravitational perturbations on the plan- 
etesimal orbits. In the bottom panel we show the azimuthally- 
averaged disk eccent ricity as a function of the distance from the 
stars, computed as in lMarzari et al.ld2.012b . for 3 different evolu- 
tionary times separated by 200 yr. The 3 curves suggest a strong 
variability of the internal disk shape with time. This is due to the 
propag ation of density wave s within the disk and was observed 
also in IMarzari et al 1 (120121) for circumstellar disks in binaries. 
The behavio ur shown in Fig. [2] lower panel, differs from that 
illustrated in IPelupessv & Zwartl (1201 3l) but the two disk mod- 
els are substantially different. IPelupessv & Zwartl (1201 3l) adopt 
a locally isothermal equation of state with a fixed temperature 
profile T(r) ~ 300 x R ~ 3/4 while our initial temperature profile 
is T(r) ~ 630 x R 1 and then it evolves in time due to radiative 
cooling and various sources of heating, including that arising 
from shock waves. In Fig. Q~|we co mpare the temperatu r e pro- 
file of the isothermal disk model by IPelupessv & Zwartl (1201 3l) 
with the averaged (over 200 yr) equilibrium profile of our radia- 
tive model. The difference in both absolute values and slope are 
noteworthy and this justifies the significant differences in terms 
of disk eccentricity. 

It has been shown in IMarzari et al.l (1201 2l) that the temper- 
ature profile and the adopted energy equation have strong ef- 
fects on the disk eccentricity in perturbed disks. A different ther- 
modynamical model may lead to significant differences in the 
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Fig. 1. Azimuthally and time averaged (200 yr) temperature pro- 
file of our radiative model at t — 10000 yr (continuous red line) 
compared with the initial t — yr non-equilibrium one (dashed 
green line) and the isothe rmal temperature profile adopted by 
IPelupessv & Zwartl d2013l) (dotted blue line). 



disk evolution and in particular in the eccentricity which may 
change by more than a factor 10. The temperature profile of 
our radiative disk i s sign ificantly higher compared to that of 
IPelupessv & Zwartl (1201 3b and this plays in fav or of a higher 
disk ec centricity in our model. On the other hand. iMarzari et afl 
dun have shown that wave propagation through adiabatic 
compressions and expansions may be more efficiently damped 
in presence of radiative losses. This might cause a reduction of 
the disk eccentricity. All these physical phenomena that can lead 
to a different disk shape are more important thatn the disk self- 
gravity which would have less influnce on the disk eccentric- 
ity than the thermodynamical model. In addition, self-gravity 
would not significantly affect the planetesimal dynamics in the 
inner regions due to the large value of the Toomre parameter Q. 
The initial different temperatur e profile is not the on l y diffe rence 
between our model and that of IPelupessv & Zwartl (120131) . also 
the superficial density is difference as theirs declines as r _I while 
ours follows a r 1 ' 2 law. 

In a real disk, not only viscous heating and radiative cooling 
determin e the disks thermal s tructure, but also stellar irradiation. 
Recently, iBitsch et"ail d2013l) have shown that stellar irradiation 
dominates in the outer regions of a disk while viscous heating 
rules close to the star where the structure of disks with and with- 
out stellar irra- diation are similar. They find via numerical mod- 
eling that only beyond approximately 8 AU from the star stel- 
lar irradiation begins to have some influence on the evolution 
of a disk. However, a significant flaring of the disk is observed 
only for massive disks with Eo ~ 3000g/cra 3 and low viscosity 
(a ~ 0.001). In this case the flaring appears consistent beyond 
30 AU from the star. In our model, the circumbinary disk of 16 
Kepler extends out to 10 AU so it is expected not to be strongly 
influenced by stellar irradiation but to be dominated by the vis- 
cous evolution. In addition, in circumbinary disks density waves 
excited by the tidal gravity field of the binary propagate within 
the disk. In three dimensions these waves act like fundamental 
modes which correspond t o large surface dist orsions in the disk 
(iLubow & Ogilvielll998l) . iBolev et all (120051) have also shown 
that shock waves in 3-D disks cause sudden increases in the disk 
scale height, a phenomenon called hydraulic jump. If the spiral 
waves excited by the binary perturbations are also shock waves, 



the hydraulic jumps in the vertical direction can shield the outer 
disc from stellar irradiation. Disk self-shadowing due to spiral 
density waves may strongly reduce the relevance of stellar irra- 
diation in circumbinary disks. 

The relatively fast changes of the disk eccentricity and of 
the disk gravity field has an additional perturbative effect on the 
planetesimals trajectories favouring eccentricity excitation. The 
orientation of the disk changes slowly with time on a timescale 
longer than 10 4 yrs. However, the computation of the planetes- 
imal orbits is very time consuming and the timespan of each 
model is limited by the amount of CPU requirement. This is also 
the reason why we neglect the self-gravity of the disk in our 
simulations. A model without planetesimals has however shown 
that self-gravity does not significantly affect the shape and time 
variability of the circumbinary disk adopted in our simulations. 
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Fig. 2. In the upper panel we show isocurves of the disk density 
distribution (in normalized units) after 2 x 10 5 binary revolutions 
illustrating the eccentric shape of the disk. In the bottom panel 
the azimuthally-averaged disk eccentricity is drawn as a function 
of the radial distance at 3 different evolutionary times separated 
by 200 yr. 



3.2. Orbital evolution of 5 km size planetesimals 

The three major sources of perturbations on the planetesimal or- 
bits are: 

- - The secular and short term perturbations of the gravity field 
of the binary 
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- - The gas drag force 

- - The disk gravitational force 

The first kind of perturbation is accounted for also in pure 
N-body models (neglecting the mutual gravitational attraction 
of planetesimals) and its effects are summarized by Eqs. 1 and 2. 
The second term is instead handled very differently in a purely 
axisymmetric approach and our model. If the axisymmetric ap- 
proximation is adopted, the gas velocity is assumed to point al- 
ways in the tangential direction respect a circle centered in the 
binary baricenter and with radius equal to the planetesimal os- 
culating radial distance. Its modulus is the local Keplerian ve- 
locity reduced by a factor that accounts for the pressure term 
v g = vjc(1 - 2rf) 1 / 2 where r/ is of the order of 10~ 3 . In our model, 
the gas velocity is computed directly from the solution of the 
hydrodynamical equations. This is an important aspect since the 
gas velocity, due to the disk elliptic shape and the presence of 
spiral waves, is very different from that computed in the axisym- 
metric approximation. The gas velocity is no longer circular and 
it can have a significant radial component. Its modulus can be 
much larger than the value estimated by the previous simplified 
formula. In addition, the cylindrical symmetry is lost and the di- 
rection and modulus of the gas velocity depend on the azimuthal 
angle. 

The third perturbative component on the motion of planetes- 
imals, that we will show to be the dominant one, is due to the 
gravity field of the disk. An asymmetric distribution of mass 
causes a significant perturbation of planetesi mal trajectories. A 
simila r phenomenon was also observed by iNelson & Gressell 
d2010l) in fully turbulent disks where embedded planetesimals 
develop large mutual encounter velocities due to stochastic grav- 
itational forces caused by turbulent density fluctuations. In our 
scenario the density has a shaped asymmetric pattern and this is 
a worse source of perturbation since its effect does not average 
to 0. In addition, the ori entation of the perihelia play an impor- 
tant r ole in our scenario (Thebault et al. 2006; M arzari & Scholll 
l2000l) . in particular when large eccentricities are excited. The 
gravitational perturbations of the eccentric disk significantly de- 
crease the level of perihelia alignment of same size planetesimals 
and this has important consequences on the accretion process. 

In Fig. [3] we illustrate the distribution of the eccentricity and 
perihelion longitude vs. semimajor axis for 5 km size (radius) 
planetesimals. A peak in eccentricity with a value around 0.4 is 
observed in between 2-3 AU with planetesimals quickly drifting 
inwards. The semimajor axis drift rate is in fact strongly depen- 
dent on the eccentri city of the body in the equation given by 
lAdachietalJ(fl976l) describing the effect of the gas drag: 



da 
dt 



77 = 



T drag 

where 

Vkep ~ 



-Or + o e + 



2° 1X( ' + i6 



e 2 + 



1 



gas 



Vkep 



(9) 



(10) 



measures of the amount by which gas orbits the star (or star 
pair) more slowly than a solid body due to the gas's partial pres- 
sure support with Vke P being the local Keplerian velocity. The 
timescale Td mg is given by: 



Tdrag — 



8pfl 



ICoPgasVkep 



(ID 



The region developing large eccentricities is then depleted 
on a short timescale. The irregular shape of the disk in the inside 
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Fig. 3. In panel 1 we show the distribution of the eccentricity 
vs. semimajor axis of 5 km size (radius) planetesimals at two 
different times. In panel 2 the distribution is that of the perihelion 
longitude vs. semimajor axis. 



region is at the origin of the large eccentricity values that lead to 
the fast inward drift. Farther out, the eccentricity is lower but still 
high compared to the pure N-body predictions. In this region the 
radial drift is reduced. The pericenter longitude is only partly 
aligned and the phasing depends on time since it changes with 
time. Only the region between 3 and 5 AU seems to maintain 
some level of coherence over time. This coherence leads to lower 
impact velocities. 

3.3. The influence of the disk gravity unveiled 

By inspecting Fig. [3] the first question that comes to mind is: is 
it the radial component of the gas drag force or the gravitational 
attraction of the disk that is responsible for the large values of the 
planetesimal eccentricities? To answer this question we run an 
additional model where the gravity of the disk on planetesimal is 
switched off. Fig.[5]shows the difference in the two cases after 1 x 
10 3 yrs of evolution. The case with the gravity of the disk acting 
on planetesimals show much larger eccentricities than the test 
case without the disk gravity. Note that the eccentricity profile in 
the model where the disk perturbations are included differ from 
that shown in Fig. [3] since the evolutionary times are different 
(10 4 yr in Fig.|3]and 10 3 yr in Fig . [5j ■ 

While it is possible to analytically estimate the effects of a 
non-linear gas drag due to an eccentric precessing gaseous disk 
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dBeauge et al.ll2010l) . it is a prohibitive task trying to predict ana- 
lytically the details of the gravitational perturbation of an asym- 
metric disk. The gravity field felt by each planetesimal is a com- 
bination of the tidal field of the binary stars, of the eccentric disk 
and of its spiral arms which are tightly wound close to the stars. 
Planetesimals are well embedded in the disk and then they are 
sensitive not only to the overall shape of the disk but also to its 
fast time variability. In Fig.|4]top panel we show azimuthally av- 
eraged radial profiles of the potential produced by the disk at dif- 
ferent evolutionary times. To give an idea of the variation of the 
potential with azimuth we plot also the variation of the potential 
with azimuth respect to the local average value at two different 
times. The inner zone, within ~ 1.8 AU, show a slow decrease 
of the potential with a limited azimuthal variability. This is the 
region where the planetesimals are less excited in eccentricity. 
Just beyond 2 AU, the potential begins to rise and the dynamics 
of planetesimals reacts to this trend change with a steep increase 
in eccentricity. It is possibly the combination of radial and az- 
imuthal variations that account for the sudden raise in the plan- 
etesimal eccentricity, even if it appears difficult to analytically 
predict the amount of perturbation felt by the planetesimal tra- 
jectories. This because the whole shape of the potential changes 
with time. In addition, there is no analytical expression avail- 
able for the gravity field of an elliptic disk within itself. While 
the outside potential might be fitted with the analytical expres- 
sion derived by McCoullogh dMurrav & Dermottlll999l) . inside 
the disk the task appears much more complex. Not only a sec- 
ular theory predicting the eccentricity perturbations of an ellip- 
tical disk on bodies located inside the disk itself is intrinsically 
very complex but to make the task even more difficult the disk 
changes with time and, as a consequence, also the potential, as 
shown in Fig. [4] 

Even if the effects of the non-radial component of gas drag 
are not relevant in exciting large eccentricity values and the disk 
gravity does all the job, this does not mean that the gas drag can 
be neglected when modeling the planetesimal evolution. Large 
values of eccentricity powers up the gas drag effect on the semi- 
major axis a since da/dt is proportional to e, as discussed in the 
previous section. As a consequence, the eccentricity excitation 
leads also to a fast inward migration related to the eccentricity 
value. Does the non-radial component of gas drag contributes 
at later times to excite the eccentricity to the values observed in 
Fig. |3]-a ? When the eccentricity is excited by the disk gravity, 
the large value of 77 is mostly due to the radial velocity induced 
by the eccentric orbit of the planetesimal rather than the irregu- 
lar value due to the the disk. As a consequence, we expect again 
that the dominant term is the disk gravity and that gas drag still 
acts as a damping force. This is further confirmed by the case of 
25 km size planetesimals discussed in the next section. It has to 
be noted that the simulations shown in Fig.[5]show the eccentric- 
ity value after 10 3 of evolution while those in Fig.|3]and Fig. [6] 
illustrate the eccentricity distribution at later times (10 4 ) when 
the planetesimal evolution reaches an almost stationary state. 

3.4. Evolution of 25 km size planetesimals 

As a further test that the disk gravity is responsible for the large 
eccentricity values of the planetesimals, we ran an additional 
simulation for larger bodies 25 km in radius. In Fig. |6]we com- 
pare the eccentricity and pericenter distribution of the two differ- 
ent size swarm. The eccentricities are much larger and in effect 
5% of the bodies are injected in hyperbolic orbits. Their orbits, 
located in between 2-3 AU where the eccentricity excitation is 
the highest, are perturbed by the disk until they reach an eccen- 
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Fig. 4. Azimuthally averaged radial profile of the disk potential 
(given in normalized units) at different evolutionary times sam- 
pled every 500 yrs. Some change in the potential is also due to 
the mass loss through the inner and outer borders of the compu- 
tational grid. The two bottom plots illustrate the azimuthal varia- 
tion, computed as AV = (V - V)/V, at t = 2500 and t = 4000 yr, 
respectively. The azimuthal variation is not regular and depends 
on the evolutionary time. 



tricity of about 0.6-0.7 so that, at pericenter, they come close to 
the binary. Repeated pericenter passages further pump up their 
eccentricity because of the interaction with the binary until a 
close encounter with the secondary star ejects them out of the 
system. In addition, the pericenters are not phased at all for large 
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Fig. 5. Test run to verify the impact of the disk gravity on the 
planetesimal dynamics. The outcome of the full model is com- 
pared with that of a model without the gravitational attraction 
of the disk on planetesimals. In this second case the forced ec- 
centricity is significantly reduced. The comparison is performed 
after 1000 yrs of evolution and for 5 km size planetesimals. 



planetesimals. This is due to the reduced effect of the gas drag 
which was partly able to damp the eccentricities of the smaller 
5 km size planetesimals but it is unable to perform this task for 
25 km size bodies. As a consequence, the disk gravity is more 
effective for the large planetesimals in exciting their orbits. 

3.5. Impact velocities 

The dynamical behaviours identified in the previous sections 
have to be interpreted in terms of how they affect mutual accre- 
tion of planetesimals. The outcome of mutual planetesimal colli- 
sions depends on both their impact velocities and their respective 
sizes. Ideally, one would thus like to follow a whole population 
of planetesimals with a given size distribution and record all col- 
lisions for all impacting pairs of sizes Ri and /?2- Unfortunately, 
because of the high CPU cost of computing planetesimal orbits 
in this sophisticated set-up, only 400 can be followed simultane- 
ously. This is not enough to consider a size distribution amongst 
them and have enough statistics on dvR u R 2 everywhere in the 
disk. We thus considered 2 simplified cases where all planetes- 
imals have the same size, one for "small" planetesimals with 
R = 5 km, and one with a "large" planetesimals with R = 25 km. 
The consequence of this simplification is to underestimate im- 
pact velocities among planetesimals, as gas drag tends to mini- 
mize these velocities for equal-sized i mpactors while increa sing 
them for differentially-sized objects dThebault et alj|2006b . As 
such, our estimates should be considered as a best case (that is, 
accretion-friendly) scenario. 

To derive an estimate of the mutual impact velocities be- 
tween the "small" and "large" planetesimals orbiting the bi- 
nary from the orbital parameters of the test bodies in our hybrid 
model, we used a post-processing code computing all possible 
crossings among the test planetesimals of our sample. For each 
pair of orbits, the code looks for the crossing location and com- 
putes the relative velocity from 2-body keplerian formulas. In 
this way we build up a statistical sample of possible impact ve- 
locities which characterize the planetesimal swarm around the 
binary. These velocities are then compared to the critical ve- 
locity v C rit(R), corresponding, for impacts between equal-sized 
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Fig. 6. In panel 1 we compare the distribution of the eccentricity 
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bodies of size R, to the limiting value above which impacts re- 
sult in net mass erosion in stead of mass accretion. To derive 
v C rit(R), we use Eqs. 1 and 2 o flLeinhardt & Steward (1201 2b . which 
lead to v cr it(5km) ~ 25 m s _1 for 5 km size planetesimals and 
v C rit(25km) ~ 125 m s _1 for 25 km size planetesimals. 

In Fig. [7] we show the relative velocities among equal size 
R = 5km (upper panel) and R = 25km (lower panel) planetesi- 
mals. We see that, for both cases, large values of impact speed 
are excited in the proximity of the inner border of the disk, i.e., 
below ~ 2 AU. This high velocity regime for equal-sized bodies 
is the direct consequence of the gas disk gravity, because pure 
gas drag would have lead to perfect orbital phasing for a given 
planetesimal size, and thus very low Av for same size bodies. 
However, the effect of gas drag is still noticeable, as it is respon- 
sible for the fact that average Av are higher for 25 km planetes- 
imals, ~ 600 - 1200 m s 1 , than for 5 km objects, ~ 200 - 400 
m s _1 . This difference can be unambiguously attributed to gas 
drag, as it is the only size-dependent mechanism acting on plan- 
etesimals, so that its damping effect on the eccentricity will be 
much more pronounced for smaller bodies. However, even for 
the small planetesimal run, velocities are still much higher than 
the critical value for erosion v cr it(5km) ~ 25 m s _1 . This is also 
true for 25 km planetesimals, despite of the fact that v C rit($km) is 
higher, ~ 125 m s _1 . 



7 



F. Marzari et al.: Planetesimal accumulation in 16 Kepler B 



1000 



100 



Q. 
E 




4 6 
Radial distance (AU) 



10 



1000 



100 



CO 
Q. 

E 




4 6 
Radial distance (AU) 



10 



Fig. 7. Relative encounter velocities between planetesimals. In 
panel 1 we compare the relative velocities in the case of small 
planetesimals at different times. There is a significant time vari- 
ability in the impact velocity which, however, remains large 
all the time. In panel 2 we illustrate the relative velocities for 
R = 25km planetesimals at t — I0 4 yr. The black lines in both 
plots show the average impact velocity computed in small radial 
bins The 2 horizontal lines display the limit ing Av values for the 
accreti on/erosion frontier, as deduced from lLeinhardt & Stewart! 
d2012l) . for 5 km impactors and 25 km ones. 



The high-Av regime is maintained, for all planetesimal sizes, 
in the strongly-depleted region between 2 and 3.5 AU. It is fol- 
lowed by a region, between ~ 3 .5 and ~ 5 .5 AU, of lower impact 
velocities for the small planetesimal case. However, these veloci- 
ties remain above the erosion threshold, except for a very narrow 
region around 4 AU, where < Av >~ v cr j,. For the larger 25 km 
planetesimals, impact speeds are > 300 m s _1 , and thus higher 
than Vent, in this whole intermediate region. 

In the outer regions of the simulations, up to 10 AU, im- 
pact velocities become comparable for both planetesimal popu- 
lations. Their values remain relatively high, in the ~ 100-200m 
s _1 range, and are above the erosion threshold everywhere, ex- 
cept for the outermost 8.5-10 AU domain, where < Av > become 
slightly lower than v„„ for 25 km objects. 



4. Discussion and Conclusions 

The discovery by KEPLER of planets in circumbinary orbits has 
reawaken the interest regarding the planetesimal accumulation 



process in such systems. The first numerical investigations of 
this issue adopted an axisymmetric approximation for the gas 
disk in which planetesimal are imbedded, implicitly assuming 
that the tidal force of the central binary does not significantly 
perturbs this disk. We show in this paper that this is not the case 
and, due to the presence of the companion star, the circumbinary 
gas disk becomes eccentric. This has profound implications for 
the accumulation process of planetesimals due to the excitation 
of their orbital eccentricity and partial destruction of the perihe- 
lia alignment. Our simulations where the disk is evolved together 
with the planetesimals indeed show that, adding the crucial ef- 
fect of the gas disk gravity greatly increases impact velocities 
amongst planetesimals in the circumbinary disc around Kepler 
16. For the two planetesimal populations we have considered, 
5 and 25km, the environment is globally very hostile to mutual 
accretion in the region within 10 AU. 

As noted earlier, we have considered a simplified case where 
all bodies have the same size, so that only a small fraction of all 
possible planetesimal encounters have been explored. However, 
in the regions where gas drag has a non-negligible effect on 
mutual impacting speeds, we expect impact velocities between 
differentially-sized bo dies to be even highe r than those amongst 
equal-sized ones (e.g. IThebault et al"1l2006l) . As a consequence, 
in these regions our runs probably display a best-case scenario, 
with a situation that would be even more accretion-hostile in a 
"real" system with a spread in planetesimal sizes. The r < 6 AU 
domain clearly corresponds to this case, the important differ- 
ence between hv^m and Av25fan being a clear indicator that gas 
drag is an important factor in imposing the impact speeds. In 
the region beyond 6 AU, impact speeds for 5km and 25km bod- 
ies are comparable, which could at first glance indicate that gas 
drag has a limited effect there. However, Fig J6] clearly shows 
that it is not the case, since these two populations have very 
different eccentricities and periastron in these outer regions, so 
that gas drag still has an important influence on collision ve- 
locities and the fact that Avs^h ~ Av'25jt m is simply a coinci- 
dence. We thus here also expect impact speeds to be higher be- 
tween differentially-sized objects than for equal-sized ones, so 
that our single-sized simulations here again probably correspond 
to the most accretion-friendly configuration possible. The fact 
that even in this best-case scenario the whole disc below ~ 9 AU 
is strongly hostile to accretion, with the possible exception of 
a narrow strip around 4AU, seems to indicate that this result 
should also hold for any planetesimal population with a spread 
in its size distributions. Furthermore, even if accretion was to be 
possible, it is important to stress that it cannot be as efficient as 
in an unperturbed case. Indeed, planetesimal accumulation will 
be slowed down, because of the increased encounter velocities 
due to perturbations which would probably switch off the fast 
runaway-growth mode that r equires very low impact velocities 
to proceed (see discussion in IThebau lt 201 1). 

Our results should, however, be taken with some caution, as 
one potentially important effect is not accounted for in our runs, 
that of the possible re-accretion of collisional fragments that 
has been identified by Paardekooper et al.(2012). This mech- 
anism would probably help planetesimal growth in this dy- 
namically excited environment, so that its omission overesti- 
mates the erosive behaviour of the disk. Unfortunately, it is im- 
possible to implement this effect in the already very complex 
set-up that has been considered here. We note, however, that 
other simplifications of our numerical approach might have the 
opposite effect, i.e., underestimating the erosive nature of the 
planetesimal swarm. The first one is that we suppose that all 
planetesimal start on circular orbits at the same time, whereas 
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there m ight be a spread in their formation epoch (e.g IXie et al] 
I2010bl) . Such a spread would lead to increased differential or- 
bits, and thus impact spee ds, even between equal-sized objects 
dPaardekooper et aI1l2012l) . The other simplification is that the 
spatial resolution in the hydrodynamical model required to cor- 
rectly compute the disk-planetesimal interactions limits the ra- 
dial extension of the disk we can model with a reasonable CPU 
load. Our gas disk is thus probably too small, and we would ex- 
pect a larger and more massive disk to have more powerful grav- 
itational perturbations on the planetesimal trajectories, possibly 
increasing the impact velocities. 

Given these uncertainties, it is thus too early to reach defini- 
tive conclusion regarding the precise balance between accretion 
and erosion in a realistic planetesimal disk. However, to the very 
least our simulations have shown that disk gravity plays a crucial 
role and always act towards increasing impact speeds and the 
erosive behaviour of the swarm. As such, they strengthen and 
expand the results obtained by previous works, which already 
identified that the region where the planet is located is hostile 
to pla netesimal accretion dMeschiarill2012l iPaardekooper et ail 
l2012h . This seems to rule out the in situ formation of the Kepler 
16 planet following the core-accretion scenario. This result also 
probably holds for the Kepler-34, Kepler-35 and Kepler-47 plan- 
ets, given the similarities between these different systems. 

A possibility to explain the present position of these plan- 
ets is that they formed farther out in the circumbinary disk. 
Subsequent migration due to the interaction with the disk would 
have brought them back to th eir present orbits, as shown in 
dPierens &Nelsodl2007ll2008h . This is also suggested by their 
mass which is in t he Neptune-Saturn range, in ag reement with 
the prediction of dPierens & Nelson! 120071 120081) . Jupiter size 
planets in fact would be either ejected from the system or sent 
on outer orbits. Our simulations show that such migration would 
have to be very efficient, bringing the planet where it is today 
from an initial formation region probably located beyond 6 AU 
from the binary's center of mass. 

An alternative scenario may be based on the direct formation 
of large planetesimals from the accumulation of small solid par- 
ticles in turbulent struct ures of the gaseous disk ([Johansen et al.l 
l2007tlCuzzi et al.ll2008l) . In circumbinary disks the onset of tur- 
bulence may be favored by the tidal gravity field of the central 
stars and this might lead to the formation of planetesimal large 
enough to sustain the high velocity impacts occurring in the in- 
ner regions of the disk. In this case, planets might form closer to 
the center of the disk by-passing the critical phase of small body 
accretion. However, these instability-based scenarios still need 
to be quantitatively investigated in dynamically perturbed envi- 
ronments such as binaries before any conclusion can be reached. 
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